%=============================================================================
%   This program computes summary statistics for the data and estimates the
%   fractional integration parameter 
%
%   ALL OF THE RESULTS ARE IN TABLE 1 IN THE PAPER
%
%   Daniela Osterrieder and Peter Schotman
%   The Volatility of Long-term Bond Returns:
%   Persistent Interest Shocks and Time-varying Risk Premiums
%
%   REVIEW OF ECONOMICS AND STATISTICS
%==============================================================================

%*******************************************
%INPUT: The program loads
% (1) rawdata.txt
% (2) ELW_ShiPhi06.m
% (3) LW_ShiPhi06.m
%OUTPUT: Matlab tables
% (a) "DATA" (table with clean data)
% (b) "TABLE1" (Table 1 in the paper)
%*****************************************

clear 
clc
close all

%importing the raw data (data runs from 1954.01 until 2012.01)
%the structure of the data file is day | month | year | T-Bill | Bond 60 | Bond 120
data=load('rawdata.txt'); 

%definition of the size
[n,col]=size(data);

%FRED FED quotes discount basis, but we want continously compound monthly yields 
t_bill=(-400*log(1-((data(:,4))/400)))/12; 
%change in t_bill
tb1=t_bill(2:end,1)-t_bill(1:end-1,1);
%Bond yields are percentages per annum on a continuous compounding basis; devide by 12 to get monthly yields
bond5=(data(:,5))/12; 
bond10=(data(:,6))/12;

%compute excess bond returns
rx_61=zeros(n,1);
rx_121=zeros(n,1);
for a=2:n
    rx_61(a,1)=-60*bond5(a,1)+61*bond5(a-1,1)-t_bill(a-1,1);
    rx_121(a,1)=-120*bond10(a,1)+121*bond10(a-1,1)-t_bill(a-1,1);
end;
clear a

%calculate spreads
s_60=bond5-t_bill;
s_120=bond10-t_bill;

%make a data table 
DATA=table(data(:,1),data(:,2),data(:,3),t_bill,[0; tb1],rx_61,rx_121,s_60,s_120,...
    'VariableNames',{'Day' 'Month' 'Year' 'TBill' 'TBill_Change' 'Excess_Return_5_Year' 'Excess_Return_10_Year' 'Spread_5_Year' 'Spread_10_Year'});

clear data n col t_bill tb1 bond5 bond10 rx_61 rx_121 s_60 s_120

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% ESTIMATION OF d  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%we estimate d for the T-Bill, and the spreads
YY=[DATA.TBill DATA.Spread_5_Year DATA.Spread_10_Year];
[T,cc]=size(YY);

%initialize
d_ELW=zeros(cc,1);
exitflag_ELW=zeros(cc,1);

d_LW=zeros(cc,1);
exitflag_LW=zeros(cc,1);

%size of the spectral window
m=round(T^0.5);

for lal=1:3
    
    x=YY(:,lal);

    %Generate the vector lambda, i.e. the harmonic (or fundamental) frequencies
    lambda=((2*pi)/T)*[1:1:T]';
    
    %%%%%%%%%%%%%% 1. Exact local Whittle %%%%%%%%%%%%%%%%%%%%%%%%%
    
    %The exact LW, which has been proven to be consistent and asy. normal for any d
    if lal==1;
        x0=x-x(1,1); %shimotsu mean correction
        d0_ELW = 0.95;
    else
        x0=x-mean(x); %shimotsu mean correction
        d0_ELW = 0.4;
    end;
    options = optimset('Display','off','TolFun',1e-15,'TolX',1e-15,'TolCon',1e-16);
    [d_ELW(lal,1),fval_ELW,exitflag_ELW(lal,1)]=fminsearch(@ELW_ShiPhi06,d0_ELW,options,m,lambda,x0,T);
    
    clear x0 d0_ELW
    
    %%%%%%%%%%%%%%% 2. Local Whittle %%%%%%%%%%%%%%%%%%%%%%%%%
    %For the stationary area, i.e. d between (-1/2 and 1/2), the local Whittle (LW) is consistent
    %and asymptotically normal.
    
    if lal==1;
        %we have to consider the TBill rate in first differences
        x0=zeros(T,1);
        for t=2:T
            x0(t,1)=x(t,1)-x(t-1,1);
        end;
        clear t
    else
        x0=x;
    end;
    
    %Now define the discrete Fouries transform of the series x (evaluated at
    %each frequencies lambda)
    w_LW_pre=(T/sqrt(2*pi*T))*exp((2/T)*1i*pi*[0:1:T-1]').*(ifft(x0));
    %change order (as in Shimotsu and Phillips, 2006)
    help=w_LW_pre(1,1);
    w_LW=[w_LW_pre(2:end,1); help];
    clear help w_LW_pre
    
    %Use the dft's to construct the periodogram of x at each frequency lambda
    I_LW=w_LW.*(conj(w_LW));
    clear w_LW
    
    d0_LW = 0;
    [qre,fval_LW,exitflag_LW(lal,1)]=fminsearch(@LW_ShiPhi06,d0_LW,options,m,lambda,I_LW);
    
    if lal==1
        d_LW(lal,1) = qre+1;
    else
        d_LW(lal,1) = qre;
    end;
    clear qre
    
end;
clear lal I_LW d0_LW cc T lambda options fval_ELW fval_LW x YY x0 

se_d=(0.25/m).^(1/2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%   SUMMARY STATISTICS    %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%we provide summary statistics for T-Bill, change in T-Bill, excess
%returns, and spreads
YY=[DATA.TBill DATA.TBill_Change DATA.Excess_Return_5_Year DATA.Excess_Return_10_Year DATA.Spread_5_Year DATA.Spread_10_Year];
[n,cc]=size(YY);

%initialize
MEAN=NaN(cc,1);
STD=NaN(cc,1);
ACF_1=NaN(cc,1);
ACF_84=NaN(cc,1);

for j=1:cc
    
    if j==2 || j==3 || j==4
        %for the change in T-Bill and excess returns the first observation
        %is zero -> delete
        vect=YY(2:end,j);
        T=n-1;
    else
        vect=YY(:,j);
        T=n;
    end;
    
    %average
    MEAN(j,1)=mean(vect);
    
    %standard deviation
    STD(j,1)=std(vect);
    
    %autocorrelation
    ACF_1(j,1)=(1/(T-1))*ones(1,T-1)*((vect(1:end-1,1)-mean(vect)).*(vect(2:end,1)-mean(vect)));
    ACF_1(j,1)=ACF_1(j,1)./var(vect);
    
    ACF_84(j,1)=(1/(T-84))*ones(1,T-84)*((vect(1:end-84,1)-mean(vect)).*(vect(85:end,1)-mean(vect)));
    ACF_84(j,1)=ACF_84(j,1)./var(vect);
    
    clear vect T
    
end;
clear j

sum_stats=[MEAN STD ACF_1];
clear MEAN STD ACF_1 n cc YY

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%   OUTPUT TABLE 1 IN THE PAPER    %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

TABLE1=table(sum_stats(:,1),sum_stats(:,2),sum_stats(:,3),[d_ELW(1,1); NaN; NaN; NaN; d_ELW(2:3,1)],[d_LW(1,1); NaN; NaN; NaN; d_LW(2:3,1)],...
    'VariableNames',{'Average' 'Std_Dev' 'First_AC' 'EW' 'LW'},...
    'RowNames',{'TBill' 'TBill_Change' 'Excess_Return_5_Year' 'Excess_Return_10_Year' 'Spread_5_Year' 'Spread_10_Year'});
clear sum_stats d_ELW d_LW


fprintf(1,'\n\n-----------------------------------------------------------------------------------------------------\n');
    fprintf(1,'-----------------------------------------------------------------------------------------------------\n');
    fprintf(1,'\n\n                      Summary statistics                          \n');
    fprintf(1,'\n\n-----------------------------------------------------------------------------------------------------\n');
    fprintf(1,'-----------------------------------------------------------------------------------------------------\n\n');
TABLE1
    fprintf(1,  '-----------------------------------------------------------------------------------------------------\n');
    fprintf(1,'The spectral window for EW and LW estimation of d has size:  \t %f       \n', m);
    fprintf(1,'The standard error for all estimates of d is: \t %f \n',se_d);
    fprintf(1,'84th autocorrelation of the T-Bill is: \t %f \n',ACF_84(1,1));
    fprintf(1,  '-----------------------------------------------------------------------------------------------------\n');
    fprintf(1,  '-----------------------------------------------------------------------------------------------------\n\n');

clear se_d m ACF_84

